A first passage problem for a bivariate diffusion process: 
numerical solution with an application to neuroscience. 
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Abstract 

We consider a bivariate diffusion process and we study the first passage time of 
one component through a boundary. We prove that its probability density is the 
unique solution of a new integral equation and we propose a numerical algorithm 
for its solution. Convergence properties of this algorithm are discussed and the 
method is applied to the study of the integrated Brownian Motion and to the in- 
tegrated Ornstein Uhlenbeck process. Finally a model of neuroscience interest is 
also discussed. 
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1. Introduction 

First passage time problems arise in a variety of applications ranging from fi- 
nance to biology, physics or psycholo gy ([|29Ll20l 12111 and examples cited therein). 
They have been largely studied (see I127ll for a review on the subject): analytical 
fl S E3, HE H , numerical or approximate results [Hffl&fifiHQSIB 
exist for specific classes of processes such as one dimensional diffusions or Gaus- 
sian processes. On the contrary, the case of bivariate processes has not been 
widely studied yet. Indeed, results are available only for specific problems such 
as the first exit time of the considered two-dimensional process from a specific 
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surface HToL [Till . However there is a set of instances where the random variable 
of interest is the first passage time of one of the components of the bivariate pro- 
cess through a constant or a time dependent boundary. Examples of this type of 
problems are the First Passage Time (FPT) of integrated processes such as the 
Integrated Brownian Motion (IBM) or the Integrated Ornstein Uhlenbeck Process 
(IOU). Indeed, these one dimensional processes should be studied as bivariate 
processes if the Markov property has to be preserved. Recent examples of ap- 
plications of the IBM or of the IOU processes have appeared in the metrological 



literature 11 1 311 where these processes are alternatively used to model the error of 
atomic clocks. In that case the crossing problem corresponds to the first attain- 
ment of an assigned value by the atomic clock error. Another application for this 
type of problems arises in neuroscience for the study of two- compartment models 



II 1511 . Indeed, the membrane potential evolution of two communicating parts of the 
neuron, the dendritic zone and the soma, can be depicted by a two-dimensional 
diffusion process, whose components describe the two considered zones. Further- 
more, the time of a spike, i.e. the time when the membrane potential changes its 
dynamics with a sudden hyperpolarization, is described as the FPT of the second 
component through a boundary. Motivated by these applications, we consider the 
FPT of one component of a bivariate diffusion process through an assigned con- 
stant boundary, we prove an integral equation for this distribution and we propose 
a numerical algorithm for its solution. In Section 2 we introduce the notations and 
the necessary mathematical background. In Section 3 we present the new integral 
equation and the condition for the existence and uniqueness of its solution. In 
Section 4 we introduce a numerical algorithm for its solution and show its con- 
vergence properties. In Section 5 we illustrate the proposed numerical method 
through a set of examples, including the two-compartment model of a neuron. Fi- 
nally in Section 6 we compare computational effort and reliability of the proposed 
numerical method with a totally simulation algorithm. 

2. Notations and Mathematical Background 

Let X(f) = (Xi(t) ,X2(t))' ', t > to, be a two-dimensional diffusion process on 
I CM 2 originated in X(?o) = y, where the superscript ' denotes the transpose of 
the vector. Let s < t, we denote with 

f(x,t\y,s) = ^A_P(X(f) <x\X(s)=y) (1) 

its transition probability density function, where x = (xi,X2) and y = (yi,^)- 
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In this paper we are concerned with the random variable FPT of the first compo- 
nent of the process X(f ) through a boundary S > y\ : 

T = mf{t>t :X l (t) >S}. 

To describe T we will use its probability density function 

g(f|y,f ) = ^P(r<f|X(fo)=y). (2) 

Furthermore let Z(t) be a random variable whose distribution coincides with the 
conditional distribution of X 2 (t ) given T = t 

F(X 2 (T)<z\T = t;X(t )=y). (3) 

We will also consider the joint distribution of (X 2 (T), T) and its probability den- 
sity function 

d 2 

g c ((S,z),t\y,t ) = —F (X 2 (T) < z, T < t \X(t ) = y ) , z e R, t e [t , «] . (4) 

We denote with Kx(h(X)) the expectation with respect to the probability measure 
induced by the random variable X. We skip the subscript if there is no possibility 
of misunderstanding. In the following theorem we link (0Q) with ©. 

Theorem 1. Forx\ > S it holds 

P(X(f)>x|X(f )=y) (5) 
rt f+°° 

= d# g c ((S,z),#\y,to)F(X(t)>x\Xi(&)=S,X 2 (&)=z)dz. 

If the joint probability density function f (x,t \ y, to) exists, it holds 

ft f+°° 

f(x,t\y,t )= / dV / g c ((S,z),&\y,t )f(x,t\(S,z),&)dz. (6) 

J tQ J— CO 

Proof. Equation © is a consequence of the strong Markov property, as explained 
in the following. 

Let h : (S, °°) x R — > R be a bounded, Borel measurable function and let J^y- be 
the c-algebra generated by the process X(f ) up to the random time T. We get 

E[h(X(t))\X(t )=y] = E[E[h(X(t))\& T ;X(t )=y}] (7) 

= E[E[h(X(t))\X(T)]] 

> r+oo 

d$ / E[h(X{t))\X(fi) = (S,z)]g c ((S,z),fi\j,to)dz 



where the first equality uses the double expectation theorem while the second one 
uses the strong Markov property. Choosing h(y) = /{ Xl! oo}x{x 2 ,°°}(y) we § et ©• 
Finally, writing the conditional probability P(X(/) > x|Xi(#) = S,X2(&) = z) as 
a double integral, changing the order of integration and differentiating © with 
respect to x\ and X2 we get ©. □ 



Remark 1. Equation © was introduced without proof in |J10f| . 

The transition probability density function (Q3 is known in a few instances. 
One of these cases is a process solution of linear (in the narrow sense) stochastic 
differential equation [0] 



dX(t) = [A(t)X(t) +M(t)]dt + G(t)dB(t), t > t 
X(f )=y 



(8) 



where A(t) and G(t) are 2 x 2 matrices, M is a vector of 2 components and B(f) 
is a bivariate standard Brownian motion. 

The solution of ([8]), corresponding to the initial value y at instant to, is 



x(0 = <Kmo) 



y+ / 0(k,/q) 1 M(u)du+ [ <j)(u,to) l G(u)dB(u 



(9) 



where (j)(t,to) is the solution of the homogeneous matrix equation 
j t <j>(t,t )=A(t)(l>(t,to) with 0(? o ,*o) = I. 



(10) 



For t E [0,o°], the diffusion process has a two-dimensional normal distribution 
with expectation vector 



m(t \j,t ):= E(X(0 |X(/ ) = y) = <^ ?o) 
and 2x2 conditional covariance matrix 



y+ / (f>(u,to) l M(u)du 
J t 



(11) 



Q(t\y,t ) = $(t,to) 



(j)(u,to) l G(u)G(u) ((j)(u,to) l ) du 



-W. 



(j)(t,to), (12) 



where the superscript ' denotes the transpose of the matrix. 
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In the autonomous case (A(f) = A, M(/) = M and G(t) = G), expressions (fTTT) 
and (PT21) are simplified 



m(t\y,t Q ) 



y+ f e-^'-'^Mdu 

Jtn 



(13) 



Q(t\y,to] 



f e -Mu-t ) GG ' e -\'(u-t ) du 
Jtn 



A (u-t ) 



(14) 



f e Mt-u) GG ' e A'(t-u) du _ 

Jtn 



For Gaussian or constant initial condition, the solution of © is itself a Gaussian 
process, frequently known as Gauss-Markov process. 

Examples of Gauss-Markov processes are the Integrated Brownian Motion (IBM), 
the Integrated Ornstein Uhlenbeck Process (IOU). Also the underlying process of 
the two-compartment neural model [|15ll is a Gauss-Markov process. 



If det Q(t | y, to) ^ for each t, the probability density function / (x, t | y, to) of any 
two-dimensional Gauss-Markov process is 



exp{-J[x-m(*|y,f )]'Q(*|y,fo) 1 [x-m(f |y,f )]} 



27ty/detQ(t\y,t ) 
and follows the Chapman-Kolmogorov equation ||2lll . 



3. An Integral Equation for the FPT Distribution 

Let us consider a diffusion process {K(t),t > 0} originated in y = at to = 0. 
It holds 

Theorem 2. If 

F(X l (t)>S\X l ($) = S,X 2 ($)=z) ,zGl,^[0,i] (15) 

and its derivative with respect to t are continuous in < # < t, then the FPT 
probability density function is the solution of the following integral equation: 

P(Xi(0>S|X(0) = 0) (16) 
= jf d# g (0 I 0, 0) E m [P (Xi (0 > S |Xi (0) = 5, X 2 (0) )] 
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where the distribution ofZ(-&) is given by ©. 
The solution of rfTol ) exists and it is unique. 

Proof. Let us consider © with x\ = S and X2 = — °°, we get 

P(Xi(0>S|X(0) = 0) (17) 

ft r+°° 

= d& gc ((S,z),#\0,0)F(X 1 (t)>S\X l (#) = S,X 2 (&) = z)dz. 
Considering 

F(X 2 {T) < Z ,T <t\X(0)=0)= f JTP(X 2 (T)<z|r = T,X(0) = 0)g(T|0,0), 

Jo 

taking the derivatives with respect to z and t, using ©, we get 

g c ((S,z),t\0,0) = -^-F(X 2 (T)<z\T = t;X(0) = 0)g(t\0,0). (18) 

Substituting (U~8l) in (fT7l) we get the integral equation ([ToT l. It is a first kind Volterra 
equation with regular kernel 

k(t,$)=E zW [F(X 1 (t)>S\X l ($)=S,X 2 ($))], 

because k(t, #) is bounded. In particular k(t, t) does not vanish for any t > 0. 
Due to the hypothesis (fT5l) . the kernel of the Volterra equation (fT6l) and its deriva- 
tive respect to ? are continuous for < # < t. 

Similarly, the left hand side of equation (TT6b and its derivative respect to t are 
continuous for t > 0. Furthermore P (Xi (0) > S |0, 0) = 0. 

Thus, applying Theorem 5.1 of (11611 . we get the existence and uniqueness of the 
solution. □ 

Corollary 3. The first passage time probability density of a Gauss-Markov pro- 
cess (HI) satisfies the following equation 

l-^l^ftl= (19) 



= jfd0g(0|O,O)E zW 



I -Erf 



S-mW(t\ (S,X 2 (fl)),fl) 
y/2QW(t\(S,X 2 (*)),*) 
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where m^\t) =m^\t\ 0. 0) denotes the first component of the vector 4771). Q^\t) 
Qt (t\ 0,0) denotes the element on the upper left corner of the matrix ([72]) and 
Erf{x) denotes the error function ^J. 
The Volterra equation rf79l) admits a unique solution if 



d f S-mM{t\{S,z),$) 

is a continuous function oft>~&>0. 
Proof. Due to the Gaussianity of the process, we have 

dx 2 J s d Xl f(x,t\y,tQ) 

= l (\ Erf( S - mil){tlyJo) 

Replacing this result into (fT6l) . we obtain ([T9i 
Since (|2Q|) is continuous for hypothesis, then 

iPWO>WA,)= y ) = ii(i-Erf(^aL)) 



(20) 



= 1 ?::r J ( S-mW(t\y, k) ) \ 2 \ d ( S-mW(t\y,t ) 

is continuous. Thus applying Theorem |2] we get the existence and uniqueness of 
the solution of (fT9l). □ 



Remark 2. The term 

F(X 1 (t)>S\Xi(#) = S,X 2 (#) = z) 

represents the probability of being over the threshold S after a time interval (t — 
$), starting from the threshold itself. For a Gauss-Markov process it becomes 



P(Xi(f) >S|Xi(0) = S,X 2 (#)=z) = \ 



1 - Erf 



S-mW(t\ (S,z),fl) 
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Thus under weak conditions of its well-definition, applying the l'Hopital's rule, 
the following limit 



*-"{ W2Q< n >(l\i>,(S,z))J J 

assume a positive value C < 2. Then using the dominated convergence theorem 
we can conclude that 

= C. (21) 

Remark 3. Note that the random variable X2(#), that appear in the expectation 
(ED, has values on an interval [a,b] that changes depending on the features of the 
process ©. 



1 - Erf 



/ S-mW(f|(<>,X 2 (fl)),fl) 
W2G(")(r|(S,X 2 (0)),0) 



4. Gauss-Markov processes: a numerical algorithm 

The complexity of equation (fl9l) does not allow to get closed form solutions 
for g. Hence we pursue our study by introducing a numerical algorithm for its 
solution. 



Let us consider the partition to — < t\ < . . . < = t of the time interval [0, t] 
with step h = — t^-i for k= l,...,N. 

Discretizing integral equation (fT9l) via Euler method, we have: 



1-ErffcfM 



(22) 



£|(0-|0,0)E z(fj) 

7=1 



1-Erf 



5-mW(^|(5,X 2 (0)),Q) 
2QM(t k \ (S,X 2 (tj)),tj) 



fork=l,...,N. 



Equation d22l) gives the following algorithm for the numerical approximation g (t |0, 0^ 
ofg(T|0,0),TG(0,0. 
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Step 1 



«e,io,o) = ^ 



1 - Erf 



where C is given by (I2TT) 
Step k, k > 1 

1 



1(^1 0,0) 



ch i J 



(23) 



1 



k-\ 



-^g(0|0,0)E z(f . ) 



1 -Erf 



S -mW(t k \(S,X 2 (tj)),tj) 
2Qt n Kt k \(SMtj)),tj) 



Note that the first term on the r.h.s. is obtained for j = k. 

To sum up, the FPT probability density function in the knots to,t\,...,tN is the 
solution of a linear system Lg = b where 



1-Erf 



and 



/ C/i 
^ d N ,ih 6 N , 2 h ■■■ 



( g(ti I 0,0) \ 

V g(t N \o,o) j 



Ch j 



with 



e kj = E z(0 .) 



1-Erf 



S-mW(t k \(S,X 2 (tj)),tj) 
2QW(t k \(S,X 2 (tj)),t/. 



(24) 
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for k = 1 , . . . , N and j = 1 , . . . , k. 

To evaluate the expected value (l24l) for k = 1 , . . . , N and j = 1, we make 
use of the following Monte Carlo method. 

We repeatedly simulate the bivariate process until the first component crosses the 
boundary and we collect the sequence {Z ; , z = 1, . . .M} of i.i.d random variables 
with probability distribution function © with t = tj. Then we compute the sample 
mean 

S-mW(t k \(S,Zi),tj) 



LEiErf 



y/2QM(t k \(S&Uj) 



e kj = i ^ ; '-. (25) 

k ' J M v 

Here M is the sample size. 

The following theorem proves that this algorithm converges. In order to sim- 
plify the notations of the theorem, let us first define 

I S-mP)(t k \ (S,z),tj) 
\jf(z;t k ,tj) := Erf im ' 

\^2Q^)(t k \{S^tj) 
for k = 1 , . . . , N, j = 1 , . . . , k. Then it holds 

Theorem 4. If the sample size M for the Monte Carlo method is such that the 
error \X\= h 2 at a confidence level a and if there exists a constant a, such that 
for all h> 

max E z( , ) [\if(X 2 (tj),t k ,tj) - \if(X 2 (tj),t k _ h tj)] < ah, (26) 

L<k<N,l<j<k— 1 

then the error £k = g (tk \ 0, 0) — g (t k \ 0, 0) of the proposed algorithm at the dis- 
cretization knots t^ for k = 1,2,... is \£k\ = 0(h) at the same confidence level 
a. 

Proof. The Euler method and the Monte Carlo method applied to (fl9l) give 

1 _ E rf ( S - mW ^) \ = y h§ t t | oo) q (27) 



while (IT9T ) can be rewritten as 

/ S-mW(t k ) 



1-Erff 5 ) = l%fe|0,0) {e k j + t,) + 5(h,t k ) (28) 
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where 8(h,t k ) denotes the error of Euler method and A indicates the error of the 
Monte Carlo method at confidence level a. 
Subtracting (1281 ) from (127T ) we obtain 

8(h,t k )=h J £(e k . J e J + Xg(t J \0,0)). (29) 

7=1 

Differencing (l29l and using (12TI) . we get 

5(/ l ,^)-5(/ l ,^_ 1 )=ft£(4 J -4-u) e ;+ /lCe * + ^(^l ' ) ;l 

7=1 

or equally 

7=1 

Then, due to the hypothesis (|26|) and to the large number law, choosing M large 
enough, we have 

1 x c-/, >. i ah k ^} g(t k \Q,0)\X\ 

\e k \ < 0-\8(h,t k ) - 5(M*-i)l + -£ £ j+ hc • 

Finally, observing that the error of Euler method is 1 8 (/z, t ) \ = 0(h 2 ) , choosing M 
such that the error of the Monte Carlo method is | A | = h 2 and applying Theorem 
7.1 of n 1611 . we get the thesis. □ 

Remark 4. In the autonomous case, hypothesis (1261) is verified. 
Indeed 

2 rP((s,z),t k ,tj) 2 
xif(z;t k ,tj)-xif(z;t k _ h tj) = -7=/. / > (30) 

< a x [p ((S,z),t k ,tj) -p ((S,z),t k - h tj)] 

where 

P {{S,z),ti,ti) - 



y/2& n Kti\{S,z),ti) 
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Using equation (fT4l) . we get that 

Q( n \t k \ (S,z),tj) = G (11) (^i I (S,z),tj) + [* e Au MM e A '"du. 

J 

Hence inequality (f30b becomes 

y/[z,tk,tj) -\\r\z,tk-\-,tj) <«2 , , (-31) 

^fi( 11 >ftt-ilOS,z),f/) 

where £?2 is a new constant. 

Using equation (PT31) . y = (5,z) and t k — t k -\ = h, k= l,...N, we can conclude 
that 

™(tk\y,tj) -m(tk-i\y,tj) 

< a 3 |(A(? A .-^_ 1 ))y + ^ A (I+A(t k -u))Mdu- ' (7+A(fjt_i - m))MJmJ 

A---2 



a 3 <j(Ay + M)/*-A(fc-j--)M/i 



Therefore inequality (1311 ) becomes 

yf(z,t k ,tj) -xif{z,tk-i,tj) < , = = 

5. Examples 

In this section we apply the algorithm presented in Sections |4] to some ex- 
amples. Firstly we consider an IBM and an IOU Process. Recent examples of 
application of the IBM and the IOU processes have appeared in the metrologi- 
cal literature [U8I1 to model the error of atomic clock. Then we will consider the 



two-compartment model of a neuron [15], whose underlying process is a bivariate 
Ornstein Uhlenbeck process. 

5.1. Integrated Brownian Motion 

The Integrated Brownian Motion by itself is not a Gauss-Markov process be- 
cause it is not a Markov process. However we can study this one dimensional 
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process as a bivariate process together with a Standard Brownian motion, as fol- 
lows 

dXi(t) = X 2 (t)dt 

(32) 

dX 2 (t) = dB t , 
with X(0) = 0. 

The process (l32l) is a particular case of the Gauss-Markov process ([8]), where 



A(t) 



1 




,M(t) 



and G(t) 




1 



There exist analytical solutions of the first passage time problem of the in- 
tegrated component of the process (|32l . but they are not efficient because they 
involve multiple integrals [13] or suppose particular symmetry properties l|loll . 
Hence, we numerically solve the first passage time problem for an IBM using the 
algorithm proposed in Section |4] A first attempt in this direction was discussed in 
Il25n . In this instance the FPT probability density function through a boundary S 
in the knots to, t\ , . . . is solution of a linear system Lg = b where 



/ 



1 - Erf 



V6S 

W 1 



\ 



and 



1 _ Frf ( ^> s 

wv / 

/ 2h 

e 2 ,ih 2h 
d3,ih 9 X2 h 2h 

\ 6 Nj ih Nj2 h ■■■ 



\g(^|0,0) / 



2h 



with 0, 



%i 



1+Erff ^ X2{ ' j) 



for k,j=l,...,N. 



Note that in this case the constant C defined in Remark |2] is equal to 2 and the 
range of the random variable X 2 (T) is [0, °°] . 
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In Figure \T\ we show the FPT probability density function of an IBM through a 
boundary 5 for three different values of the boundary. 



.... S=1 

■-. S=3 

? | S=6 

0.25 - 



0.2 - ; : 

o : ■ 

g 0.15 - ; ■ 

o) : : 



: I I 

o.i -; i i 




5 10 15 20 25 30 35 40 



Figure 1 : Evaluation of the FPT probability density function for an IBM through 
three different boundaries: 5=1 (dotted), 5 = 3 (dashdot), 5 = 6 (solid). 

5.2. Integrated Ornstein Uhlenbeck Process 

As the IBM, the IOU Process is not a Markov process and it should be studied 
as a bivariate process together with an Ornstein Uhlenbeck Process, as follows 

dXi(t) = X 2 (t)dt 

(33) 

dX 2 (t) = (-aX 2 (t)+n)dt + adB t , 
with X(0) = 0. The process (1331) is a Gauss-Markov process d8), where 
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Note that in this case the constant C defined in Remark |2] is equal to 1 and the 
range of the random variable Xi(T) is [— °°, °°] . 

In Figure [2] we show the FPT probability density function of the IOU through a 
boundary S = 6, for /i = 0.01, o = 1 and three different values of the parameter 
a. 




Figure 2: Evaluation of the FPT probability density function of the IOU through 
a boundary 5 = 6 for /! = 0.01, o = 1 and three different values of the parameter 
a: a = 0.01 (dotted), a = 0.3 (dashdot), and a = 0.5 (solid). 



Note that the curve for a = 0.01 in Figure |2]is very similar to the curve for S 
in Figured] indeed if /i — > and a — V IOU converges to a standard IBM. 

5.3. Two-compartment model 



One dimensional neuronal models [|2l|l identify the membrane potential val- 
ues on the different parts of the neuron with those assumed in the trigger zone. To 
improve the model (11511 proposes a two dimensional approach. The membrane po- 
tential of a neuron is described by means of a bivariate stochastic process, whose 
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components represent the depolarizations of two distinct parts, the trigger zone 
and the dentritic one. 



Let (-X'i(0)f>o an d (X 2 (t)) t >o ^ e me stochastic processes associated to the depo- 
larization of the trigger zone and the dendritic one, respectively. Then, assuming 
that external inputs, of intensity /i and variability o, influence only the second 
compartment and taking the interconnection between the parts into account, we 
obtain the following stochastic model 

dX 1 (t) = {-aX 1 (t) + a r [X 2 (t)-X 1 (t)]}dt 

(34) 

dX 2 (t) = {-aX 2 (t) + a r [Xi (t) -X 2 (t)]+n}dt + odB t 

with X(0) = and where a and a r are constant related to the spontaneous mem- 
brane potential decay and to the intensity of the connection between the two com- 
partments, respectively (cf. Fig [3]>. 




Figure 3: Schematic representation of the two-compartment approach 

The process (l34l) is an Ornstein Uhlenbeck two-dimensional process, particular 
case of the Gauss-Markov process ([8]) where 

f-a-a, a,- \ = 1 0\ = f \ 

w \ a r -a-a r J w \ H J V 0cr / 

Note that in this case the constant C defined in Remark |2] is equal to 2 and the 
range of the random variable X 2 (T) is [fcS,°°], where 

a + a r 



Assuming that after each spike the system is reset to its initial value, the time be- 
tween two spikes, i.e. the time when the membrane potential changes its dynamics 
with a sudden hyperpolarization, is described by the FPT of the first component 



16 



through a boundary S. 

In n 1 5M this model was studied using simulation techniques, but applying the algo- 
rithm proposed in Section |4] we can compute the interspike intervals distribution 
as the FPT probability density function of the first component of X(f). In Figure 
@]we show the FPT density for a set of values of the input. 
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Figure 4: Evaluation of the FPT probability density function through a boundary 
5 = 6 for a two-compartment stochastic model of a neuron, choosing a = 0.33, 
a r = 0.2, o = 1 and three different values of fi: jl = \ .5 (solid), /i = 2 (dashdot), 
and }i = 3 (dotted). 



6. Numerical algorithm versus a totally simulation algorithm 

The introduced numerical method involves a Monte Carlo estimation to eval- 
uate the expected value (l24l) . One may wonder about the advantages of the pro- 
posed method with respect to a totally simulation algorithm. Indeed it is easy to 
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simulate M sample paths of the considered bivariate process to get a sample of 
FPTs. These FPTs could be used to draw histograms or their continuous approxi- 
mations. However this approach is computationally expensive. Indeed it requires 
large samples to give reliable results. Moreover the estimation of the tails of the 
distribution is scarcely reliable and time consuming. 

On the contrary the numerical method proposed here requests weak compu- 
tational efforts despite the presence of the Monte Carlo method and is applied to 
estimate a specific integral, i.e. the expectation (l24l) . Indeed, also a small number 
of trajectories guarantees reliable results. 

In Figure |5] we compare the accuracy of the results obtained with the two 
methods. We simulate M sample paths of the IBM in order to determine a sample 
of M FPTs and we use it to draw the corresponding histogram. The same sample is 
used to compute the sample mean (1251) to get the FPT probability density function 
via the numerical method. The choice M = 1000 gives reliable results in both 
cases. However, when M = 50 the histogram is rude while the numerical method 
does not loose its reliability. We further underline how the computational time 
to build the histogram or to draw the FPT density with the proposed numerical 
method are comparable, for the same value of M. 

In Figure |6] we show the shapes of the FPT probability density function ob- 
tained via the proposed numerical algorithm. We use a sample of size M = 50 
(solid line) and M = 1000 (dash line) to compute (|25l ). Their differences are neg- 
ligible. 

7. Conclusion 

We studied the FPT of one component of a bivariate diffusion process through 
a boundary. We wrote a new integral equation for the FPT probability density 
function proving its existence and uniqueness. A numerical algorithm for its solu- 
tion was developed proving its convergence properties. The algorithm was applied 
to a set of processes of interest for various applications. Advantages of the method 
with respect to a totally simulation algorithm are discussed. 

The crossing problem for one component of a multivariate process or the case 
of random initial value can be treated as extensions of the proposed equations and 
of the numerical algorithm. 
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Figure 5: FPT probability density function for the IBM obtained via the proposed 
numerical method and the corresponding histogram. Sample of size M has been 
used to build the histogram and to evaluate (|25l) in the numerical method: (a) 
M = 50 (b) M = 1000. 
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Figure 6: FPT probability density function for the IBM obtained via the proposed 
numerical method computed using a sample of size M = 50 (solid line) and M = 
1000 (dash line) for the sample mean (|25l ). 
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